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Abstract. This paper describes an algorithm for determining radii of convergence of power expansions for 
algebraic functions and the testing done to check it. Since the current methods for computing these series 
are iterative, standard methods for computing radii of convergence cannot in general, be used. However, 
relying on geometric properties of algebraic functions, convergence radii of these series can be determined 
precisely. 



1. Introduction 

The objects studied in this paper are fractional power expansions of algebraic functions 

f(z, w) = a (z) + dx(z)w + a 2 (z)w 2 H + a n (z)w n = 0, (1) 



i — isuch that f(z,w) is irreducible over the rationals with z and w complex variables and the coefficients, 
cii(z), polynomials in z with rational coefficients. The degree of the function is the highest power of w. 

. By the Implicit Function Theorem, this equation defines locally, an analytic function w(z) when — — ^ 0. 

^ ® w 

+^ And by Newton- Puiseux's Theorem[7], Q can be factored over the field of fractional power series, K ({z}) 

f(z,w) = f[(w-w d (z) j ), w d (z) G K ({z}) , (2) 

Q\ where each Wd{z) is a branch, possibly multivalued, of the function. These fractional power series have 
^ a radius of convergence at least equal to the distance to the nearest singular point. However, one result 
O °f this paper is to demonstrate the actual radius of convergence can and often does extend beyond many 
l/-^ singular points. In the test cases described below, one branch is shown to extend across 118 singular points. 
O Since the current methods for computing fractional power series of algebraic functions rely on iterative 
^ means [8J, radii of convergence, in general, cannot be calculated by the ordinary techniques such as the 
r| ratio test or root test. This paper describes a method to compute these values precisely. 
• The purpose of this paper is four-fold: 

(1) Describe in detail, the Newton polygon algorithm as described by Kung pE], 

(2) Demonstrate that a numeric version of Newton Polygon can be implemented with satisfactory 
results, 

(3) Implement an algorithm for determining the radius of convergence of an algebraic power series, 

(4) Provide a software tool for further investigating algebraic functions. 

2. Some properties of algebraic functions used in this paper 

Fractional power expansions of algebraic functions are called Puiseux series and are in this paper com- 
puted by the method of Newton polygon. These series often represents only a small portion of the function 
near the origin. However, the entire function can be represented by Puiseux series in annular discs sur- 
rounding the origin. These annular expansions can in principle, be computed by numeric means using a 
modified version of Laurent's expansion theorem and if computed accurately, can be used in the algorithms 
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below to determine their annular domain of convergence. However, this paper focuses only on power series 
inside the disc D(0, \r c \). 

The resultant of f(z,w) with it's derivative f w is denoted by R(f,f w ). A point where one or both of 
a n (z) or R(f, f w ) is zero is a singular point of /. A point where a n (z) = is also a pole, possibly ramified, 
of the function. In this paper, singular points are labeled as s„. A property of algebraic functions which 
distinguishes them from single-valued functions is algebraic functions can be both singular and analytic at 
a singular point. That is, a singular point may not affect all coverings of an algebraic function unless the 
function is fully-ramified at the singular point. For example, a 10-degree function may have only a single 
2-cycle branch and eight 1-cycle coverings at a singular point. In this case, the 2-cycle covering is singular. 
The eight single-cycle coverings are not analytically affected at this singular point unless one is affected 
by a pole of the function. However, if the function fully-ramifies into a 10-cycle branch at this singular 
point, all coverings would be affected. It is for this reason algebraic power series can have a radius of 
convergence extending beyond the first singular point: Their branch coverings may simply not be singular 
at a singular point. Only when the covering becomes singular does the convergence radius of its power 
expansion become established. The main objective of this paper is to identify which singular point is 
interrupting the analyticity of branch cycles thereby establishing the radius of convergence of their power 
expansions. 

The following conventions are used in this paper: 

(1) The Puiseux expansions of the roots of consist of a set of <i-valued branches. A branch is 
sometimes called a d-cycle where d is a positive integer. A power series in z 1//3 would be a 3-cycle 
branch. It has three coverings over the complex z-plane just like the function f(z) = ^fz. The sum 
of the cycles is always equal to the degree of the function in w. 

(2) In order to identify a particular branch, w(z) of /, a two-level identification is used. The first 
level is the cycle number representing the number of coverings. The second identifier is the sort 
order of the branch: A base point z m on the real axis is selected mid-way between the origin and a 
disc of radius equal to the distance to the nearest non-zero singular point. When the function has 
multiple d-cycle branches, the values {u^(z m )} are computed and ordered by increasing real part 
then imaginary part. The first value in the sort order for each branch is identified and the branches 
are then ordered according to this sort order as 1^,1,^,2 an d so forth. 

(3) The concept of "branch" is used throughout this paper and differs in meaning to the definition 
encountered in basic Complex Analysis. In this paper, the term "branch" refers to a multi-valued 
d-cycle root of f(z,w) given by (j2l). 

(4) The following discussion makes use of the term, "extending a branch over a singular point". This is 
in reference to the discussion above about singularities and the coverings they affect. 

(5) In this paper, r c is a positive integer representing a singular ring when non-zero singular points are 
arranged in order of increasing absolute value. The smallest non-zero singular point is therefore 
on ring one, singular points with the next largest absolute value are on ring two and so forth. |r c | 
however, represents the absolute value of the associated singular points. 



3. Newton polygon method 

Given ([T|, we seek to compute power expansions 

.. 00 

= iE^T ( 3 ) 

A/ 

n=0 

where E is the normal exponent and d is the cycle number, for all branches of the function. The simplest 
method of computing ^ is by the method of Newton polygon. For a function of degree n in w, the method 
computes a set of n series representing the roots of However, the actual branches of the function are 
often multivalued and ramify or wind around the origin multiple times. For example, a ten-degree function 
can ramify into a single-valued branch, a 2-cycle branch, a 3-cycle, and 4-cycle branch. Four of the power 
series computed by the method would be conjugates of the 4-cycle, three series, conjugates of the 3-cycle, 
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FIGURE 1. Flowchart for Newton Polygon 



and so forth. Thus in the case of the 4-cycle branch, the algorithm produces the set 

W^} = ^E c 4e W4 )"(- 1/4 ) n , J =0,1, 2, 3. (4) 

n=0 

In the particular implementation of the algorithm presented here, a basis set of power series is produced 
consisting of a single power series for each o?-cycle branch. The remaining conjugate series can be computed 
from the basis series using the appropriate form of Q. 

A theoretical basis for the Newton polygon algorithm can be found in Walker [5J. The particular 
implementation of the algorithm for this paper is taken from Kung [4J, a flowchart of which is shown in 
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Figure [T} We begin first by representing a branch as 

w(z) = c x z Xl + c 2 z Xl+X2 + c 3 z Xl+X2+X[i + ■■■ (5) 
in which Ai > and Aj > for i > 1 and q 7^ 0. The steps in the process are as follows: 

(1) Normalizing the function: 

In order to achieve the requirement above on the exponents Aj, the function is normalized to 
eliminate negative exponents in the expansion. The normalization process removes algebraic poles 
at the origin from the function which arise when zero is a root of a n (z). For example, 

f(z, w) = (78 + A9z - 46z 2 ) + (-2z)w + (91 - 3Az - 80z 2 )w 2 + (52z + A7z 2 )w 3 

has a pole at the origin and thus one branch of the function would have a negative exponent in its 
expansion. The normalization process produces for this function, 

g(z,w) = z 2 f(z,w/z) 

= (78z 2 + A9z 3 - 46z 4 ) + (-2z 2 )w + (91 - 3Az - 80z 2 )w 2 + (52 + A7z)w 3 , 

which in this case E — 1. We can then expand g(z, w) around the origin, which now has no singular 
point there and hence no expansions with negative exponents, and then multiply each series by the 
quantity 1/z to arrive at the expansions for f(z,w). 

(2) Computing the Newton polygon and extracting the lower Newton leg: 

The principle device for beginning the computation of the coefficients c n for the series is the 
Newton polygon which for the function f(z,w) is the convex hull of it's support. The support of 
the function is a set of points, {i,a(cii)} where a(oj) is the order of cii(z). The order of cii(z) is 
the lowest power of z in ai(z) with the exception that if ai{z) = 0, then a(cii) = 00 and the point 
(i, 00) is omitted. Drawing the convex hull of these points and taking the left-most lowest segments 
so that all remaining points are either above or to the right of the segments then gives the lower 
Newton leg. The Newton polygon for 

f(z, w) = (z 14 + 3z 15 ) + (2z 15 + 2z 16 )w + (3z 15 - 20z m )w 2 

+ {Az 15 )w 3 + (10z 5 -z 6 + 2z 7 )w A + {8z w )w 5 + {9z 10 )w 6 + {20z)w 7 (6) 

+ (3z)w 8 + (2)w 9 + (5)w w 

is shown in Figure [2j In this example, the lower Newton leg has four segments given by the red 
lines. The blue dashed lines provide the slope and vertical intercept for each segment which gives 
the exponents (3 and A discussed below. 

(3) Deriving the characteristic equation: 

Choose a segment from the lower Newton leg and form the characteristic equation of the segment. 
The principle is based on letting w(z) = z Xl (c± + wi) and then making the substitution f(z, z Al (ci + 
Wi)). Extracting the lowest power of z alone and setting it to zero then provides a means of 
computing c\. However, the form of the characteristic equation can be read off directly from the 
polygon diagram: for each segment, choose the points intersected by the segment. In the case of 
Figure [2j the top most segment intersects the points (0, 14) and (4, 5). These are from the order-14 
term of do(z) and the order-5 term of a^(z). Now choose the coefficient of each of those terms and 
label them &degree,order = ^>i,a- I n the two cases above, we have lo,i4 and 10^5. We can write this set 
of points as: 

Now form the characteristic equation 

Q 

K{x) = Y, b w xij =° ( 7 ) 

with Q being the total number of points. This give for the segment above, K(x) = 1 + 10x 4 . The k 
distinct, non-zero solutions to the characteristic equation are the first coefficients C\ of k power series 
for the function. If the zeros are multiple roots, then we create a series for each set of multiple roots 
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Figure 2. Newton polygon for m 



by setting c\ for each root, and let f\(z, w±) 



-ft 



1 f(z, z x (ci +Wi)) where Ai is the negative of the 



slope for the segment and /3i, the vertical intercept of the segment. We then compute a second-level 
polygon for f\(z,Wi) with the exception that segments with zero slopes are now omitted from the 
Newton leg, derive a second characteristic equation for C2, and continue in this way if necessary for 
C3 and so forth until we reach fk{z,Wk) with simple roots for its characteristic equation. 
(4) Convert to regular form: 

Once simple roots are obtained for the characteristic equation, we can use a variation of Newton 
iteration to compute additional terms of the series. In order to use this method, we need a poly- 
nomial with integer powers; the procedure above will often produce an expression with fractional 
powers. But we can convert the expression fk(z,Wk) above to one with integer powers with the 
following two substitutions: 

f(z,w) = -^-f k (z,z Xk w) 
z* h 

J(z,w) = f(z d ,w) 



(8) 



where k is the last recursion of the function which did not produce a characteristic equation with 
multiple roots, and d is the lowest common denominator of the exponents {Xj}, j = 1, 2, • • • , k. 
(5) Normal Iteration: 

The terms in f(z,w) not containing a factor of z make up the characteristic equation in w so 
that the first simple root computed above, r^, is a zero to f(0,w). That is, /(0, r^) = 0. Thus, we 
let in the iteration step below, po — r fc- We can see this with f(z, w) = (97z 7 ) + (100z 2 + 77z 3 — 
77z A + Q4z 5 )w + (94 — 7bz 2 )w 2 + (— 54)w 3 . The first segment produces a characteristic equation of 
97 + lOOx. When we process the function as above, f(z, w) = 97 + 100u> + zg(z, w) and therefore, 
the root of the characteristic equation is a root of f(0,w). The iteration process is a variation of 
Newton iteration, a description of which can be found in Kung. We let 



Pn+l =Pn~ mod —=— 

V df(p. 

po = n, 



z 2 " +1 



(9) 
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where the mod function is simply extracting the 2™ +1 — 1 terms of the Taylor series of the quotient. 
After a set number of iterations we obtain the n-th iteration and write for the branch, 

i 

3=1 
t 

e* = X>, ( 10 ) 
k=i 

1 ' 

1=1 

where E is the normal exponent of the function. In this way, the method produces n power series 
for an n-th degree function. 

4. Implementing a numerical version of Newton Polygon in Mathematica 

Existing implementations of Newton polygon use exact arithmetic due to the possibility of numerical 
errors that produce incorrect polygon results. However, these exact methods limit the types of functions 
that can be studied due to the computational complexity, storage space and execution times involved. In 
this paper, a numeric version of the method is implemented and shown to produce acceptable results for 
a variety of different function types. Numerical approximations however introduce the following issues: 

(1) Multiple root resolution, 

(2) Residual coefficients close to zero, 

(3) Loss of accuracy due to coefficient expansion or contraction, 

(4) Numerical integration drift. 

These issues are dealt with as follows: 

4.1. Multiple root resolution and residual coefficients. Given the function: 

f(z, w) = (w- l)(w - 2f{w - 3) 3 (u> - 4) 4 - z, 
when we create the polygon for this equation, we obtain as the characteristic equation for c±, 

27648 - 110592a; + 192384a; 2 - 192832a; 3 + 123852a; 4 - 53428a; 5 

+ 15715a; 6 - 3118a; 7 + 400a; 8 - 30a; 9 + a; 10 = 0. (11) 

Solving this for x using machine-precision, the roots are 

{1., 2., 2., 2.9993, 3.00009, 3.00103, 3.9948, 3.99966, 4.00375, 4.00589} 

so that if we used a tolerance of just 10~ 3 , we would not pick up the multiple roots. In order to successfully 
resolve the multiple roots, the algorithm must check the roots against the accuracy of the computation. If 
the roots do not differ by this accuracy, then the algorithm used in this study considers the roots multiples. 
Mathematica has features to monitor and control the accuracy and precision of all computations. We 
therefore begin the algorithm by computing the singular points of the function with a high degree of 
precision. This then determines the initial precision of succeeding calculations. In this study, the precision 
of the singular points, {s n }, was set between 400 and 5000 digits. We then carefully monitor the accuracy 
of the computations. If the number of significant digits or precision of a computation drops to zero, then 
the quantity is considered numerically zero. Consider the two numbers: 

nl = 2.3333333333333333333345353 

n2 = 2.3333333333333333333379889. 

These numbers differ in the 21 'th decimal place. Therefore, if the accuracy or number of significant digits 
to the right of the decimal place of the numbers is set to a value greater than 21, then the precision of 
the difference nl — n2 will be greater than zero. However, if the accuracy of the numbers are both set to 
20, then the precision of the difference nl — n2 drops to zero. This is the principle used to check both for 
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multiple roots and residual coefficients in the computations. If the roots of the characteristic equation are 
computed to 150 digits of accuracy and then that accuracy is manually decreased to 145, then multiple 
roots ri and r 2 when subtracted, will produce a number with zero significant figures or zero precision. 
Likewise, if during the Newton polygon phase of the computation, the precision of the iterates, f^(z,w), 
drops to zero when the accuracy is manually decreased to an accuracy below its current accuracy, then that 
value is considered numerically zero and dropped. This of course raises the possibility of numerical error 
if significant results are below this threshold. However the algorithms designed in this paper were written 
with diagnostic tools to detect this and other possible errors due to numerical precision. One tool is a table 
reporting the difference |n — r2 1 for all roots of the characteristic equation. The user can visually inspect 
these values if there is a concern for numerical error. Another tool reports all residual terms removed from 
the computations. These can also be checked and if necessary, the precision of the computations increased. 

4.2. Accuracy reduction. One limiting factor affecting the accuracy of the results is the computation of 
roots to the characteristic equation. The accuracy of the roots is a function of the precision of the input 



data. Another factor is the quotient _ - " in the iteration phase of the calculation. After each iteration, 

df(Pn) 

the accuracy of the results is usually less than the input data. After several iterations, the accuracy can 
drop significantly. In this implementation of the algorithm, the accuracy of all results are continuously 
checked and if it goes below zero, the algorithm stops and prompts the user to increase the precision of 
the calculations. 

For complicated functions such as the twenty-degree function studied in Test |5j the Newton iteration 
process produces extremely large coefficients on the order of several hundred. Therefore, if the precision 
of the initial data was only 200 digits, and the coefficients blow-up to values exceeding 10 200 , then the 
number of significant digits to the right of the decimal place drop to zero or even negative. The solution 
to this problem is simply to increase the precision of the input data. In order to detect when this occurs, 
the accuracy of the computations are checked and when the accuracy drops below zero, the user is notified 
to increase the working precision. 

4.3. Numerical integration limitations. The function studied in Test [5] has a third-order pole located 
at the origin. During the numerical integration phase of the continuation algorithm, the starting point for 
the integration is very close to the origin. During the testing, the numerical results varied greatly from 
the expected results. Although this problem may have been addressed by careful tuning of the numerical 
integrator, another remedy is to use the normalized function for the continuation analysis as was done in 
the test since the pole at the origin is removed. Another potential problem is when the integration path 
encounters widely- varying function values or if the path is very long. In Test [6j the integration path is of 
the order of 10 26 . If a numerical problem is suspected, the precision of the numerical integration can be 
increased in an effort to resolve it. Another solution is to carefully monitor the integration as was done for 
Test [6] as described in the Conclusion section below. 

5. Determining the radius of convergence of an algebraic power series 

The geometry of algebraic functions described above is used to determine the radii of convergence of 
their associated power series. The process is depicted in Figures [3] and |4} In general, the procedure checks 
each branch of the function for analytic continuity over singular points. Figure [4] shows the setup over the 
first singular point in red which has a nearest neighbor in blue. The first step is to compute the expansion 
around the origin . This expansion will converge at least to the nearest singular point with distance Rl 
in the figure. The value of the function is then computed at one-half the distance of Rl at the point z s 
using the expansions computed above. These values are then adjusted to the actual values as determined 
by the function Nearest(/(^ s , w) = 0). The algorithm first determines the smallest difference between the 

roots of f(z s ,w) = 0. This value is given by p m i n . All comparisons are then made to a tolerance of t-^p. 

where N is an integer which in the test cases below, was 100. If this tolerance is not met at the current 
number of terms of the series, the the user is prompted to generate series with more terms. The next step 
is to numerically extend each branch from z s in a radial direction to z e . This is done using the differential 
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equation for the function with initial value at z s . These values are adjusted to the exact values at this 
point using the method described above. The third step is to compute the local expansion around the 
selected singular point. We know these expansions will converge in a radius at least equal to the distance 
to the nearest singular point or gRl. Taking one-half the distance, we compute the value of each local 
branch at z e and adjust for exact values. If the analytically-continued values from z s to z e impinge upon 
ramified branches or poles at s n , and if the singular points were processed in sequential order, the power 
series has a radius of convergence equal to this singular ring. If a branch at the origin does not impinge 
onto local singular branches for all singular points in a ring, the power series for this branch has a radius 
of convergence beyond this ring. We continue in this fashion until no further continuations are possible. 
The details of each step are as follows: 

(1) Computing the singular points {s n }: 

In this study, the singular points were computed to sufficient precision to mimimize the potential 
for possible numerical errors described above. This precision was usually between 400 and 5000. 

(2) Computing the function basis Bf. 

The function basis is the set of conjugately-distinct power series centered at the origin. Thus 
for a 4-cycle, the basis would have only one series representing this branch. Likewise for the other 
multi-valued branches. As many as 4035 terms of a series were generated in this study. 

(3) Computing Bf(z s ) and E(w s ): 

A line is drawn from the origin to the selected singular point. On this line and one-half the 
distance to the first singular point, the point z s is chosen. Using the basis computed above, we 
compute the values of the function branches at z s . These values will differ from the exact values 
(to numerical precision) of f(z s ,w) = 0, depending on how many terms of the series are used. The 
values are adjusted by comparing the computed value Wd(z s ) to the ideal values and choosing the 
closest match and then checking that the computed values are less than p min /100. The exact values 
are E(w s ). 

(4) Continue E(w s ) to the point z e : 

For each singular point, the distance to it's nearest neighbor is computed. We choose a point 
which is one-half this distance and on the line to the origin and mark a point z e . If this point is 
smaller in absolute value to z s , it is adjusted (made smaller) so that \z e \ — \z s \ > 0. Now using the 
differential equation for the function, 



we numerically continue each branch from z s to the point z e , and as was done with z s , adjust the 
values to the closest match, E(w e ). 
(5) Compute basis B g at selected singular point: 

We next compute the basis around the chosen singular point s n by making the substitution 
g(z,w) = f(z + s n ,w) and thus an expansion around zero for g(z,w) is an expansion of f(z,w) 
around s n . The singular points are chosen in sequential order according to their radial distance 
from the origin in order to avoid integrating over a singularity. This could occur if two singular 
points are in alignment with the numerical integration path and analyticity over the closest singular 
point was not checked prior to checking the more distant point. The basis, B g , will contain at least 
one singular branch. This singularity will either be a pole and or a ramified branch with multiple 
values. Poles and ramified branches are barriers to extending the radius of convergence of power 
expansions across singular points and therefore if a branch from the basis Bf impinges upon any 
of the singular branches of B g at s n , its power series has a ring of convergence equal to the ring 
number of the current ring, r n . If a branch does not impinge upon the singular branches B g for 
all singular points of a ring, then the series has a radius of convergence which extends across this 
singular point. We therefore select the singular branches from B g , compute B g (z e ) for each sheet 
of the branch, adjust those values to their exact values, and then compare those values to each of 
the continued values Bf(z e ). Any branches in which Bf(z e ) = B g (z e ) to numerical precision has 
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FIGURE 3. Branch Continuation Flowchart 



a power series with a radius of convergence equal to \s n \. Branches which do not agree are then 
checked against the next singular point. 
(6) Continue the above process until no remaining series can be extended over the next singular point. 
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FIGURE 4. Branch Continuation Path 



6. Algorithm Testing 

In order to obtain some empirical measure of the accuracy of the test results, all power series were 
checked with values z c and such that |r„_i| < \z c \ < \r n \ and |r„| < \zj\ < l r n+i| with adjustments made 
if r c is the first or last ring. If the computed ring of convergence is correct and the series are precise, the 
partial sums will exhibit convergent behavior in the former case and divergent behavior in the later case. 
Additionally, the series were checked against 100 random points in the range j^j < \zi\ < ^|^ c | with the 
values compared against the expected value given by the function Nearest(/(2;j, w) = 0). The maximum 
error of the set is reported in the tables. 

6.1. Test Cases. 

(1) Function with cycles one through four: 

f(z, w) = (z 14 + 3z 15 ) + (2z 15 + 2z w )w + (3z 15 - 20z 16 )w 2 + (4z 15 )w 3 

+ (Wz 5 -z 6 + 2z 7 )w A + (8z 10 )w 5 + (9z 10 )w 6 (12) 
+ (20z)w 7 + (3z)w 8 + (2)w 9 + (5)w 10 . 

Continuation results for function (12) are shown in Table [T] The test was run with a precision of 
800 using 54 terms of each series. Note the continuation column. In the first row, the four and 
three cycle branches are continuable over the first singular point. These series then have a radius 
of convergence at least equal to the absolute value of the second singular point or approximately 
0.3329. Note however, only the 3-cycle branch is continuable across the singular points in the 
second and third rings but not the fourth. This means 1^3 has a radius of convergence equal to the 
size of the fourth ring or approximately 0.636. 

Convergence results are tabulated in Table [2] with r c representing the ring of convergence as 
described above. The partial sum analysis was run for each branch and agreed with the computed 
radii of convergences in the table. Divergent behavior for W3 and W4 was not observed until both 
partial sums had accumulated more than 2000 terms. 
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Table 2. Radius of Convergence Results for Test [T] 
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io- 7 




1 


0.00247 


254 


10~ 6 




4 


0.6363 


4094 


10 -14 


W4 


2 


0.3329 


4086 


io- 19 



Table 3. Radius of Convergence Results for Test [2] 



Cycle 


r c 


\ r c\ 


Terms 


Max error 




61 


1.093 


1024 


io- 15 




1 


0.1168 


1024 


io- 12 


w 3 


1 


0.1168 


1024 


io- 10 




4 


0.505 


998 


io- 7 


w 5 


14 


0.6413 


4035 


io- 10 



Function having cycles 1, 2, 3, 4, and 5: 

f{z, w) = (z 30 + z 32 ) + (z u + z 20 )w 5 + (z 5 + z 9 )w 9 + (z + z 3 )w 12 + (6)w u + (2 + z 2 )w l \ (13) 

Results of this test are in table [3j This function has a 1-cycle branch which extends across 61 
singular rings containing 118 singular points. Cycles one through three were tested with 1024 
terms at 1000 digits of precision. In order to obtain definitive convergence/divergence behavior of 
the partial sum plots for the 5-cycle branch, 4035 terms were used because the partial sums initially 
exhibited convergent behavior but began diverging after approximately 2500 terms. 

Multiple n-cycle branches: 

f 2 (z,w) = zw 4 [a(l-a 2 w 2 )] 4 - [a(a 2 -w 2 )Y , a = 3 - 141 (14) 

This function has three 4-cycle branches which cannot be extended past the first singular point. 
This implies that the first singular point must have at least three ramified sheets and in fact the 
function ramifies into two 2-cycle branches (with additional single-cycle branches) at this singular 
point. Convergence results are in Table [4j Branch w 4i3 is an example of a ramified pole. Partial 
sum testing confirmed the computed rings of convergence. 

The maximum error observed for these series is relatively large. This can be explained due to the 
fact that the radius of convergence is extremely small, on the order of 10 -16 . The random points 
used to check the series were therefore all very close to the limit of convergence for each series and 
thus a relatively large difference between the series value and actual value would be expected for 
the small number of terms tested. 
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Table 4. Radius of Convergence Results for Test [3] 



Cycle 


r c 


Vc\ 


Terms 


Max error 




1 


2.96 x 10" 


-16 


257 


10~ 4 


W 4 ,2 


1 


2.96 x 10" 


-16 


257 


10" 4 


W4,3 


1 


2.96 x 10" 


-16 


128 


10" 4 



Table 5. Radius of Convergence Results for Test [4] 



Cycle 


r c 




Terms 


Max error 




1 


0.9585 


1026 


10 -io 


W&,2 


1 


0.958 


1026 


io- 12 



(4) Function with a recursive Newton polygon: 



f(z, w) = ((w 3 + z 2 ) 2 + z 3 w 2 ) 2 + z V 



(15) 



This function has two 6-cycle branches and when analyzed with Newton polygon, produces a 3-level 
recursive polygon tree which means it generates multiple roots for two polygon phases. This checks 
the algorithm's ability to detect multiple roots at a level other than the first when the coefficients 
are exact with infinite precision. Results of this test are shown in Table [5] using 1026 terms of each 
series. 
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(5) 20-degree function with third-order pole at origin 



f{z,w) 



+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 



-Uz - Q8z 4 + 83z 5 + 88z b ) 
-20 + 192 - 19z 2 + 42z 6 + 25/ + 54/ + ffiz 9 )w 
-16z + 59z 2 + 54z 3 + 36z 5 - 91z 6 - Uz 7 )w 2 
47 - 172 5 + 642 6 + 94z 8 + 68z 9 )w 3 
-62z + 272 2 - 25z 3 + 39z 4 )w 4 
-4 + 22 + 11/ - 13/ 
532 + 5/ - 65/ + 57/ - 75z 1Ch 



85/°W 5 



)w 



-13 + 922 + 232 2 + 2 
852 - 392 4 + 782 5 + 48/ 



152 5 + 232 7 )w 7 
26/ + 22 1 V 8 



2z 2 - 5Qz d 



9z 10 )w 9 



-17 / - 65 z s + 772 4 + 642 s - 452 9 + 9Qz 



w 



(16) 



39 + 84/ + 902 4 - 62 5 + 6z b - 57/ + 39/W 



ii 



12 



.10^13 



-482 + 792 2 - 222 5 + 752 6 - 32 9 )w 

55 + IOO2 2 - 582 6 - 192 7 + 832 9 - Alz 
67 + 852 - 82 s + 162 7 )w; 14 
-232 5 - 392 7 - 202 10 )w 15 

56 + 762 + 572 2 + IOO2 3 - 402 4 + 682 s - 552 7 + 502 9 + 52 : 
69 - 422 + 532 2 - 892 4 - 132 7 + 552 8 )w 17 
-122 2 - 222 8 - ll2 10 )w 18 

62 + 232 - 42 2 - 992 4 + 92 5 - 992 6 + 572 9 - 902 10 )w 19 
922 3 - 9l2 8 + 632 9 )w 20 . 



lo^ w l6 



This function exemplifies what can happen when numerically integrating near a pole. In step 4 of 
the continuation process, we numerically integrate from z s to z e . However, if z s is near a pole the 
numerical integration may suffer, as in this case, due to the large derivative involved. Although 
it may be possible to minimize this problem with efficient use of the numerical integrator, we can 
avoid the problem by continuation over the normalized function. Recall, the normalization process 
removes the singular point at the origin. When normalized though, the powers on 2 are usually 
increased. In this particular function, the largest power of 2 rises to 63. This causes the coefficients 
of the iteration functions, p n , to grow extremely large and thus, to maintain accuracy of the data, 
the working precision must be increased. In this particular case, it was increased to 1000. Table 
[7] summarizes the results of this test. The partial sum study for this function agreed with all 
convergence results. 



(6) Function of degree 55: 

10 

f(z,w) = z-H(w-jy. (17) 

i=i 

This function represents an inverse of a 55-degree polynomial and was analyzed with 128 terms at 
a working precision of 400. Also, the working precision and accuracy of the numerical integration 
was set to 40 and 30. Results of this test case are shown in Table [7] and agreed with the partial 
sum testing. 
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TABLE 6. Convergence results for Test [5] 





7" 

' c 


nr 

V c\ 


Tprrn c 
J. CI lllo 


l\/lfl"V rTTTlT 


^1,11 


Q 
O 






1 


w l,2 


o 
z 


D 1 ^QD71 

U. IOjU i 1 




m~ 8 


w l,3 


o 
z 






1U 




4 
^± 


n 9^1 




1U 


"a, 5 


A 

■4: 








"a, 6 


o 


fl 1 ^QD71 




1U 


"a, 7 


1 7 








w l,8 


1 7 






10- 7 

1 u 


in- 

"a, 9 


i 
i 






1U 


w l,10 


i 
i 






1U 


^1,12 


i 
i 






1U 


w l,13 


i 
i 






1U 


^1,14 


u 


D 9fi71 98 




m~ 8 


^1,15 




D 9fi71 98 




10~ 8 


w l,16 




D 9^7Q1 Q 




1 n~ 7 

-LU 


"^1,17 


5 


257919 


256 


10~ 6 


^1,18 


5 


0.257919 


256 


10- 7 


Wl,19 


6 


0.267128 


256 


10- 12 


™1,20 


6 


0.267128 


256 


10- 7 


Wl,l 


13 


0.530749 


256 


10" 5 


Table 7. 


Convergence Results for Test [6] 


Cycle 


r c 


l r d 


Terms 


Max Error 


Wi 


9 


2.39 x 10 38 


256 


10- 12 




8 


2.9 x 10 32 


256 


10~ 6 


w 3 


7 


4.79 x 10 26 


256 


10- 7 


W4 


6 


1.95 x 10 21 


256 


10^ 6 


W 5 


5 


3.75 x 10 16 


256 


10~ 6 


Wq 


3 


7.24 x 10 12 


256 


10" 6 


Wj 


2 


3.69 x 10 10 


256 


10~ 6 


w 8 


1 


2.12 x 10 10 


256 


10~ 6 


Wg 


1 


2.12 x 10 10 


256 


nr 5 




4 


2.34 x 10 13 


256 


10- 4 



(7) Function with continuations across single-sheet poles: 

f(z, w) = (3 + Az) + (-6z 2 - %-)w + ( \ - lQz + ^)w 2 + {--2z + 12z 4 )^ 3 

2 2 8 4 

(Z 22/2 . rj 1 Z Z Z /I. O 

15 + — + )w 7 + + + 2z 4 )w 8 . 
3 15 ; V 1000 25 2 5 ; 



18) 



The continuation table of this function is given in Table [8j The poles are in red and since all 
the branches are single-sheets, the notation Wd, n is dropped and only the sort orders are reported. 
All but u> 18 are continuable across the first two poles and branch w\ t i continues over all poles of 
the function. Convergence results are in Table [9j Cycle w 18 at 256 terms did not produce a low 
maximum error at 256 terms. This may have been due to the fact that this branch sheet is 15000 
at the origin and quickly reaches 140, 000 at the boundary of the testing data. However, at 2048 
terms, the maximum error was reduced to 10~ 9 . 



PUISEUX SERIES 
TABLE 8. Continuation results for Test [7] 



15 



I -\ m or 

rung 




c ■ 

1** 1 


( ,/"iT"i"t" i mi p r l r^n q 


i 
i 




u.uiyyuo 


11 9 ^ ZL ^ fi 71 
jl, Z, 3, 4, 0, 0, / J- 


9 


n 1 


fi 1 


11 9 ^ 4 ^ fi 71 
jl, z, o, 4, 0, 0, I] 


Q 
O 


n 9qpi n 1 ^4 




11 ^ ^ fi 71 
jl, o, 0, 0, / J- 




n 9Qfi -u n 1 ^7 




ii ^ ^ fi 7l 


4 


_0 353 - o 346z 


0.494 


11 6i 




-0.353 + 0.346i 


0.494 


{1} 


5 


0.1-0.5* 


0.5004 


{1} 




0.1 + 0.5i 


0.5004 


{1} 


6 


-0.0733 - 0.543i 


0.548 


{} 



Table 9. Convergence Results for Test [7] 



Cycle 


r c 




r c \ 


Terms 


Max Error 


W\,\ 


6 


0.5489 


256 


io- 7 


™1,2 


3 


0.3288 


256 


io- 9 


Wl,3 


3 


0.3288 


256 


io- 8 


w 1A 


3 


0.3288 


256 


io- 7 




3 


0.3288 


256 


io- 7 


Wlfi 


4 


0.4943 


256 


io- 7 


Wl,7 


4 


0.4943 


256 


io- 7 


Wl,8 


1 


0.01997 


2048 


IO" 9 



Table 10. Run-time Statistics (in seconds) 



Test Case 


Pw 


Terms 




Timing 




S 


N P 


A c 


1 


800 


54 


1.5 


0.17 


21.5 


2 


800 


17 


11.7 


0.17 


4605 


3 


2500 


32 


0.1 


40.2 


21.4 


4 


800 


66 


0.1 


9.3 


6.5 


*5 


1000 


63 


150.2 


0.8 


2871 


6 


400 


128 


0.5 


793 


4181 


7 


1500 


256 


0.8 


63 


155 



* Results for normalized function 



7. Timing Statistics 

The algorithms were run on a dual-Pentium processor machine at 2.2 GHz under Mathematica version 
8.0. Table 10 summarizes typical timing results for the test cases. S is the timing to compute singular 
points and some auxiliary tables, Np, the timing for the Newton polygon algorithm, and Ac, the timing 
for the continuation algorithm. 
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Table 11. NDSolve results for Test [6j singular point 7 



C!vHe 


^ n ppf 
O lice \j 


Ar*lnsil \/siln£> 

/1L L U.CX1 VCL1 LLC 


1\T T ) tpqiiIIq 

1 i LJ kjvjl VC ICoUllo 


III n^y^n t 


W\ 


1 
1 


1 nnnnnn 

l.UUUUUU 


i nnnnnn 

l.UUUUUU 


U. A 1U 


z 


1 


2 000000 - 00005z 


2 000000 - 00005i 


8 79 x 10 -20 




2 


2.000000 + 0.00005i 


2.000000 + 0.00005i 


8.79 x 10^ 20 


W 3 


1 


2.9637 - 0.0433i 


2.9637 - 0.0433z 


1.137 x 10~ 19 




2 


3.1006 


3.1006 


1.555 x HT 19 




3 


2.9637 + 0.0433i 


2.9637 + 0.0433i 


1.137 x 10~ 19 



8. Conclusions 

The limiting factor controlling radii of convergence of algebraic power series is branch-sheet continuation 
to either a ramified covering or single-sheet pole. This was the principle factor used in the continuation 
algorithm described in this paper. For each test case studied, the partial-sum testing agreed with the 
computed radii of convergence for all branches: the partial sums exhibited convergent behavior for |r c _!| < 
|z c | < \r c \, and showed divergent behavior in for \r c \ < \zd\ < |^ c +i|- And since the singular points are 
the roots to the resultant polynomial with rational coefficients, these values can be computed to arbitrary 
precision for a wide class of functions. 

However, the software may require tuning for a particular function. For example, the algorithms rely 
on determining an exact match for function values when these values are only computed approximately. 
The matching is done by finding the closest match to within a tolerance of ^ of the smallest separation 
between function values. The code will detect when this tolerance is not met and will prompt the user 
to either generate more terms of the series in the Newton polygon phase or to numerically integrate with 
a higher working precision. Other potential problems could come from extremely close polynomial roots, 
very small or large coefficients encountered during the analysis, algebraic branches which have their sheets 
extremely close to one another, and other errors associated with numerical integration of a differential 
equation as was shown in Test [5j Steps to identify these errors were described above and diagnostic tools 
were designed to identify possible problems. 

Consider Test [6j The singular points are extremely large and therefore, numerical integration from the 
point z s to z e has to traverse increasingly large distances. One branch extends out to 10 38 . How can we be 
certain the integration does not over-track the branch and land on a different sheet along the integration 
path? One way to minimize this is to integrate over a well-behaved section of the function. This function 
has no poles and so is well-behaved along the integration paths. Another way is to integrate with a high 
degree of accuracy. The test case used an accuracy goal of 30 digits. And still a third way of handling this 
possible problem is with diagnostic tools. The code written for this work has an option for reporting the 



results of the integration. Table [TT] is an example of one such diagnostic tool. It was generated during the 
analysis of Test [6] about the seventh singular point near —4.97 x 10 26 . The values are however reported 
to six decimal places to avoid clutter and so some numeric quantities for this particular function appear 
to either be identical or zero but are actually different at a greater precision. The important point is 
the last column which reports the difference between the actual value of the function at z e and the value 
determined by numerical integration. The table reports that after integrating a distance of approximately 
10 16 , the integration did not vary by more than 10~ 19 and since the maximum distance between the sheets 
at z e was approximately unity, this small difference compared to the sheet separations lends credence to 
the accuracy of the results. In this case, the sheets of the function were basically flat. However, had the 
sheets values varied greatly in the integration intervals, the analysis may have encountered problems. 

However, as stated earlier, one objective of this work was to produce a first version of a software tool 
that could successfully compute Puiseux series and their radii of convergence for a wide variety of functions 
not accessible by exact arithmetic means. As best as could be determined, the algorithm produced good 
results for the test cases that were studied. If all phases of the analysis are carefully monitored with the 
available diagnostic tools, the potential for errors can be reduced. Still though, the code has much room for 
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improvement and can provide a means of computing radii of convergence precisely for the types of functions 
studied in this paper. A greater variety of function types such as those with more polygon recursions, 
functions with higher ramified singular points, more complicated coefficients, or other morphologies to 
stress-test the software would produce a more robust algorithm as would a more careful analysis of the 
accuracy of the computations. 

One might consider using the partial sum testing to conceivably compute these convergence radii. Simply 
compute a highly accurate series with many terms and check it's partial sum behavior in successive rings 
until it begins to diverge. However, the point at which the series begins to diverge is not known and as 
was observed in the test, the partial sum may even appear to slowly converge and then start to diverge. 
This is what happened for the four and three-cycle branches in Test [T] 

In principle, we should be able to dispense with the numerical integration step of the continuation 
algorithm since the convergence domain of the power series being checked, and convergence domain of 
the power series centered at the chosen singular point will always overlap using the method described in 
this paper. We would then need only compare the values from two different power series. Doing this 
would eliminate the numerical integration step in the algorithm thus avoiding potential numerical errors. 
However, this would entail the possibility of computing the value of a power series at a point very close 
to it's radius of convergence and near a singularity so might require the need for many terms to achieve a 
desired accuracy. 
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